wd)options(scipen=999)
options(ggrepel.max.overlaps = Inf)
The spreadsheet includes Alison’s original sample names; we can use
this information to associate the new sample names, which are made up of
DESeq2 GLM variables, with the old names, which reflect
Alison’s wet-lab, library-prep, etc. work
featureCounts tablesfeatureCounts tablefeatureCounts tablet_fc’s positional information in a
GRanges objectpos_info will be used in DESeq2 processing,
post-processing, etc.
dds—i.e., a “master”
model matrixdds stands for “DESeq2 dataset” and is a
DESeqDataSet objectdds are
strainstatetimekit (tcn for “Tecan”, ovn
for “Ovation”)transcription (N for “nascent”,
SS for “steady state”)auxintimecoursereplicatetechnical
DESeq2DESeqDataSet, ddscounts_data for the featureCount
talliescol_data for setting up the GLMpos_info for adding feature metadata, subsequent
subsetting, etc.dds <- DESeq2::DESeqDataSetFromMatrix(
countData = counts_data,
colData = col_data,
design = ~ strain, # Vary on strain: n3-d vs o-d
rowRanges = pos_info
)
# Make a back-up of the DESeqDataSet object
bak.dds <- dds
# How do things look?
# dds %>% BiocGenerics::counts() %>% head()
# dds@rowRanges
# dds@design
# dds@assays
dds#TODO Let’s keep this in mind and try it if we come to
think lowly expressed genes are skewing results.
# threshold <- 1000
# dds_filt <- dds[rowSums(BiocGenerics::counts(dds)) >= threshold, ]
#
# Breakdown
# 0 13166
# 1 12822
# 2 12719
# 5 12540
# 10 12358
# 20 12144
# 50 11764
# 100 11403
# 200 10927
# 500 10015
# 1000 8822
#
# rm(threshold, dds_filt)
norm_non <- dds[dds@rowRanges$genome == "S_cerevisiae", ] %>%
SummarizedExperiment::assay() %>%
as.data.frame()
norm_non$feature_init <- dds@rowRanges$feature_init[
dds@rowRanges$genome == "S_cerevisiae"
]
# Associate non-normalized values with feature metadata
norm_non <- dplyr::full_join(
norm_non,
t_fc[(t_fc$genome == "S_cerevisiae"), 1:9],
by = "feature_init"
) %>%
dplyr::as_tibble()
norm_r <- DESeq2::rlog(
dds[dds@rowRanges$genome == "S_cerevisiae", ],
blind = FALSE
) %>%
SummarizedExperiment::assay() %>%
as.data.frame()
norm_r$feature_init <- dds@rowRanges$feature_init[
dds@rowRanges$genome == "S_cerevisiae"
]
# Associate normalized values with feature metadata
norm_r <- dplyr::full_join(
norm_r,
t_fc[t_fc$genome == "S_cerevisiae", 1:9],
by = "feature_init"
) %>%
dplyr::as_tibble()
More details on this relatively new method of normalization, which combines inter- and intra-sample normalization methods and (appears to) perform quite well:
# Isolate raw counts for samples of interest
raw <- dds %>%
SummarizedExperiment::assay() %>%
as.data.frame()
raw$feature_init <- dds@rowRanges$feature_init
# Associate non-normalized values with feature metadata
raw <- dplyr::full_join(
raw,
t_fc[, c(seq(1,9))],
by = "feature_init"
) %>%
dplyr::as_tibble()
# Calculate counts per kb of gene length (i.e., correct counts for gene
#+ length); gene length is initially in bp and converted to kb
rpk <- ((raw[, 1:5] * 10^3) / raw$length)
rpk[, 6:14] <- raw[, 6:14]
# Calculate normalization factors using the raw spike-in (K. lactis) counts
norm_KL <- edgeR::calcNormFactors(
raw[(rpk$genome == "K_lactis"), ][, 1:5]
)
# Create factor for categories (groups)
model_variables <- stringr::str_split(colnames(rpk[, 1:5]), "_") %>%
as.data.frame(
row.names = c(
"sample", "stage", "day", "kit", "tx", "aux", "tc", "rep", "tech"
),
col.names = paste0("s", c(1:5))
) %>%
t() %>%
tibble::as_tibble()
group <- factor(
# Second level is numerator, first level is denominator
model_variables$sample,
levels = c("o-d", "n3-d")
)
rm(model_variables)
# Create edgeR DGEList object composed of S. cerevisiae counts per kb gene
#+ length
dgel <- edgeR::DGEList(
counts = rpk[rpk$genome == "S_cerevisiae", ][, 1:5],
group = group
)
# In the DGEList object, include the normalization factors calculated from
#+ spike-in information
dgel$samples$norm.factors <- norm_KL
# Check that the normalization factors for each library are appropriately
#+ assigned
dgel$samples
# Scale the values to counts-per-million
norm_g <- edgeR::cpm(dgel) %>% tibble::as_tibble()
norm_g[, 6:14] <- rpk[rpk$genome == "S_cerevisiae", 6:14]
norm_g
# Clean up unneeded variables
rm(raw, rpk, norm_KL, group)
rm(dgel) #TODO Delete dgel? Or use it for trying out DE analyses with edgeR?
# Isolate raw counts for samples of interest
raw <- dds %>%
SummarizedExperiment::assay() %>%
as.data.frame()
raw$feature_init <- dds@rowRanges$feature_init
# Associate non-normalized values with feature metadata
raw <- dplyr::full_join(
raw,
t_fc[, 1:9],
by = "feature_init"
) %>%
dplyr::as_tibble()
# Calculate counts per kb of gene length (i.e., correct counts for gene
#+ length or do an "RPK normalization"); then, divide RPK-normalized elements
#+ by the sum of sample RPK divided by one million: this does the actual TPM
#+ normalization
rpk <- tpm <- ((raw[, 1:5] * 10^3) / raw$length)
for (i in 1:ncol(rpk)) {
tpm[, i] <- (rpk[, i] / sum(rpk[, i] / 1e6))
}
tpm[, 6:14] <- raw[, 6:14]
norm_t <- tpm[tpm$genome == "S_cerevisiae", ]
rm(raw, rpk, tpm)
# Make the following code generic --------------------------------------------
#+ ...so that we can try it with different normalization objects (counts
#+ normalized in different ways)
# norm <- norm_non
norm <- norm_r
# norm <- norm_g
# norm <- norm_t
# Create a PCAtools "pca" S4 object for the normalized counts ----------------
#+ Assign unique row names too
obj_pca <- PCAtools::pca(
norm[, c(1:5)],
metadata = dds[dds@rowRanges$genome != "K_lactis", ]@colData
)
rownames(obj_pca$loadings) <- make.names(
dds[dds@rowRanges$genome != "K_lactis", ]@rowRanges$feature,
unique = TRUE
)
# Determine "significant" PCs with Horn's parallel analysis ------------------
#+ See Horn, 1965
horn <- PCAtools::parallelPCA(mat = sapply(norm[, c(1:5)], as.double))
Warning: more singular values/vectors requested than availableWarning: more singular values/vectors requested than availableWarning: more singular values/vectors requested than availableWarning: more singular values/vectors requested than availableWarning: more singular values/vectors requested than availableWarning: more singular values/vectors requested than availableWarning: more singular values/vectors requested than availableWarning: more singular values/vectors requested than availableWarning: more singular values/vectors requested than availableWarning: more singular values/vectors requested than availableWarning: more singular values/vectors requested than availableWarning: more singular values/vectors requested than availableWarning: more singular values/vectors requested than availableWarning: more singular values/vectors requested than availableWarning: more singular values/vectors requested than availableWarning: more singular values/vectors requested than availableWarning: more singular values/vectors requested than availableWarning: more singular values/vectors requested than availableWarning: more singular values/vectors requested than availableWarning: more singular values/vectors requested than availableWarning: more singular values/vectors requested than availableWarning: more singular values/vectors requested than availableWarning: more singular values/vectors requested than availableWarning: more singular values/vectors requested than availableWarning: more singular values/vectors requested than availableWarning: more singular values/vectors requested than availableWarning: more singular values/vectors requested than availableWarning: more singular values/vectors requested than availableWarning: more singular values/vectors requested than availableWarning: more singular values/vectors requested than availableWarning: more singular values/vectors requested than availableWarning: more singular values/vectors requested than availableWarning: more singular values/vectors requested than availableWarning: more singular values/vectors requested than availableWarning: more singular values/vectors requested than availableWarning: more singular values/vectors requested than availableWarning: more singular values/vectors requested than availableWarning: more singular values/vectors requested than availableWarning: more singular values/vectors requested than availableWarning: more singular values/vectors requested than availableWarning: more singular values/vectors requested than availableWarning: more singular values/vectors requested than availableWarning: more singular values/vectors requested than availableWarning: more singular values/vectors requested than availableWarning: more singular values/vectors requested than availableWarning: more singular values/vectors requested than availableWarning: more singular values/vectors requested than availableWarning: more singular values/vectors requested than availableWarning: more singular values/vectors requested than availableWarning: more singular values/vectors requested than availableWarning: more singular values/vectors requested than available
# Determine "significant" principle components with the elbow method ---------
#+ See Buja and Eyuboglu, 1992
elbow <- PCAtools::findElbowPoint(obj_pca$variance)
# Evaluate cumulative proportion of explained variance with a scree plot -----
scree <- draw_scree_plot(obj_pca, horn = horn$n, elbow = elbow)
scree
# save_title <- paste0("panel-plot", ".", "scree", ".pdf")
# ggplot2::ggsave(paste0(args$directory_out, "/", save_title), scree)
#TODO Work up some logic for outfile names, location(s) for outfiles, etc.
# Save component loading vectors in their own data frame ---------------------
loadings <- as.data.frame(obj_pca$loadings)
# Evaluate the component loading vectors for the number of significant PCs
#+ identified via the elbow method plus two
PCs <- paste0("PC", 1:(as.numeric(elbow) + 2))
top_loadings_all <- lapply(
PCs, get_top_loadings, x = loadings, z = "all", a = TRUE
)
top_loadings_pos <- lapply(
PCs, get_top_loadings, x = loadings, z = "pos", a = TRUE
)
top_loadings_neg <- lapply(
PCs, get_top_loadings, x = loadings, z = "neg", a = TRUE
)
names(top_loadings_all) <-
names(top_loadings_pos) <-
names(top_loadings_neg) <-
PCs
# rm(PCs)
# top_loadings_all$PC1 %>% head(n = 20)
# top_loadings_pos$PC1 %>% head(n = 20)
# top_loadings_neg$PC1 %>% head(n = 20)
# Analyze positive, negative loadings on axes of biplots ---------------------
#+ Look at the top 15 per axis
images <- list()
mat <- combn(PCs, 2)
for(i in 1:ncol(mat)) {
# i <- 1
j <- mat[, i]
PC_x <- x_label <- j[1]
PC_y <- y_label <- j[2]
images[[paste0("PCAtools.", PC_x, ".v.", PC_y)]] <- plot_biplot(
pca = obj_pca,
PC_x = PC_x,
PC_y = PC_y,
loadings_show = FALSE,
loadings_n = 0,
meta_color = "strain",
meta_shape = "replicate",
x_min = -100,
x_max = 100,
y_min = -100,
y_max = 100
)
# x and y ranges for non-normalized counts
#+ x_min = -200000,
#+ x_max = 200000,
#+ y_min = -200000,
#+ y_max = 200000
#+
#+ x and y ranges for rlog-normalized counts
#+ x_min = -50,
#+ x_max = 50,
#+ y_min = -50,
#+ y_max = 50
#+
#+ x and y ranges for GeTMM-normalized counts
#+ x_min = -50000,
#+ x_max = 50000,
#+ y_min = -50000,
#+ y_max = 50000
#+
#+ x and y ranges for TPM-normalized counts
#+ x_min = -10000,
#+ x_max = 10000,
#+ y_min = -10000,
#+ y_max = 10000
#DECISION For now, go with rlog and TPM; continue to test all four
images[[paste0("KA.", PC_x, ".v.", PC_y)]] <-
plot_pos_neg_loadings_each_axis(
df_all = loadings,
df_pos = top_loadings_pos,
df_neg = top_loadings_neg,
PC_x = PC_x,
PC_y = PC_y,
row_start = 1,
row_end = 15, # 30
x_min = -0.1,
x_max = 0.1,
y_min = -0.1,
y_max = 0.1,
x_nudge = 0.02,
y_nudge = 0.04,
x_label = x_label,
y_label = y_label,
col_line_pos = "black",
col_line_neg = "red",
col_seg_pos = "grey",
col_seg_neg = "grey"
)
# x and y ranges (etc.) for non-normalized counts (messy)
#+ x_min = -0.75,
#+ x_max = 0.75,
#+ y_min = -0.33,
#+ y_max = 0.33,
#+ x_nudge = 0.02,
#+ y_nudge = 0.04,
#+
# x and y ranges (etc.) for rlog-normalized counts (nice and clean)
#+ x_min = -0.1,
#+ x_max = 0.1,
#+ y_min = -0.1,
#+ y_max = 0.1,
#+ x_nudge = 0.02,
#+ y_nudge = 0.04,
#+
#+ x and y ranges (etc.) for GeTMM-normalized counts (a bit messy)
#+ x_min = -0.5,
#+ x_max = 0.5,
#+ y_min = -0.5,
#+ y_max = 0.5,
#+ x_nudge = 0.04,
#+ y_nudge = 0.02,
#+
#+ x and y ranges (etc.) for TPM-normalized counts (a bit less messy)
#+ x_min = -0.5,
#+ x_max = 0.5,
#+ y_min = -0.5,
#+ y_max = 0.5,
#+ x_nudge = 0.04,
#+ y_nudge = 0.02,
images[[paste0("KA.", PC_x, ".v.", PC_y)]]
}
# How do things look?
images$PCAtools.PC1.v.PC2
images$KA.PC1.v.PC2$PC_x_pos
images$KA.PC1.v.PC2$PC_x_neg
images$KA.PC1.v.PC2$PC_y_pos
images$KA.PC1.v.PC2$PC_y_neg
# images$PCAtools.PC1.v.PC3
# images$KA.PC1.v.PC3
# images$PCAtools.PC1.v.PC4
# images$KA.PC1.v.PC4
# images$PCAtools.PC2.v.PC3
# images$KA.PC2.v.PC3
# for(i in 1:length(names(images))) {
# # i <- 2
# vector_names <- names(images) %>% stringr::str_split("\\.")
#
# if(vector_names[[i]][1] == "PCAtools") {
# save_title <- paste0("panel-plot", ".", names(images)[i], ".pdf")
# ggplot2::ggsave(
# paste0(args$directory_out, "/", save_title), images[[i]]
# )
# } else if(vector_names[[i]][1] == "KA") {
# save_title <- paste0(
# "panel-plot", ".", names(images)[i], ".1-x-positive.pdf"
# )
# ggplot2::ggsave(
# paste0(args$directory_out, "/", save_title), images[[i]][[1]]
# )
#
# save_title <- paste0(
# "panel-plot", ".", names(images)[i], ".2-y-positive.pdf"
# )
# ggplot2::ggsave(
# paste0(args$directory_out, "/", save_title), images[[i]][[2]]
# )
#
# save_title <- paste0(
# "panel-plot", ".", names(images)[i], ".3-x-negative.pdf"
# )
# ggplot2::ggsave(
# paste0(args$directory_out, "/", save_title), images[[i]][[3]]
# )
#
# save_title <- paste0(
# "panel-plot", ".", names(images)[i], ".4-y-negative.pdf"
# )
# ggplot2::ggsave(
# paste0(args$directory_out, "/", save_title), images[[i]][[4]]
# )
# }
# }
#TODO Work up some logic for outfile names, location(s) for outfiles
# Plot the top features on an axis of component loading range ----------------
#+ ...to visualize the top variables (features) that drive variance among
#+ principal components of interest
p_loadings <- PCAtools::plotloadings(
obj_pca,
components = getComponents(obj_pca, 1),
# components = getComponents(obj_pca, 1:5),
rangeRetain = 0.05,
absolute = FALSE,
col = c("#785EF075", "#FFFFFF75", "#FE610075"),
title = "Loadings plot",
subtitle = "Top 5% of variables (i.e., features)",
# shapeSizeRange = c(4, 16),
borderColour = "#000000",
borderWidth = 0.2,
gridlines.major = TRUE,
gridlines.minor = TRUE,
axisLabSize = 10,
labSize = 3, # label_size
drawConnectors = TRUE,
widthConnectors = 0.2,
typeConnectors = 'closed',
colConnectors = 'black'
) +
# ggplot2::coord_flip() +
theme_slick_no_legend
p_loadings
#TODO Work up some logic for saving the plot
# Evaluate correlations between PCs and model variables ----------------------
#+ Answer, "What is driving biologically significant variance in our data?"
PC_cor <- PCAtools::eigencorplot(
obj_pca,
components = PCAtools::getComponents(obj_pca, 1:5),
metavars = c("strain", "replicate"),
# metavars = c("strain", "replicate", "sample_name"),
col = c("#785EF075", "#648FFF75", "#FFFFFF75", "#FFB00075", "#FE610075"),
scale = FALSE,
corFUN = "pearson",
corMultipleTestCorrection = "BH",
plotRsquared = TRUE,
colFrame = "#FFFFFF",
main = bquote(Pearson ~ r^2 ~ correlates),
# main = "PC Pearson r-squared correlates",
fontMain = 1,
titleX = "Principal components",
fontTitleX = 1,
fontLabX = 1,
titleY = "Model variables",
rotTitleY = 90,
fontTitleY = 1,
fontLabY = 1
)
Warning: strain is not numeric - please check the source data as non-numeric variables will be coerced to numericWarning: replicate is not numeric - please check the source data as non-numeric variables will be coerced to numeric
PC_cor
# Get lists of top loadings for GO analyses ----------------------------------
# for(i in c("PC1", "PC2", "PC3", "PC4")) {
for(i in c("PC1", "PC2")) {
# i <- "PC1"
# Positive
loadings_pos_PC <- rownames(top_loadings_pos[[i]])[1:500]
save_title_pos_PC <- paste0(
"top-500.",
stringr::str_replace_all(get_name_of_var(loadings_pos_PC), "_", "-"),
".", i, ".txt"
)
# readr::write_tsv(
# dplyr::as_tibble(loadings_pos_PC),
# paste0(args$directory_out, "/", save_title_pos_PC),
# col_names = FALSE
# )
#TODO Work up some logic for location(s) for outfiles
# Negative
loadings_neg_PC <- rownames(top_loadings_neg[[i]])[1:500]
save_title_neg_PC <- paste0(
"top-500.",
stringr::str_replace_all(get_name_of_var(loadings_neg_PC), "_", "-"),
".", i, ".txt"
)
# readr::write_tsv(
# dplyr::as_tibble(loadings_neg_PC),
# paste0(args$directory_out, "/", save_title_neg_PC),
# col_names = FALSE
# )
#TODO Work up some logic for location(s) for outfiles
}
Here, we subset out the K. lactis features. Thus, we are using only S. cerevisiae features in the size-factor estimation. No control genes are used.
dds_SC <- BiocGenerics::estimateSizeFactors(
dds[dds@rowRanges$genome != "K_lactis", ]
)
dds_SC@colData
DataFrame with 5 rows and 11 columns
strain state time kit transcription auxin
<factor> <factor> <factor> <factor> <factor> <factor>
n3-d_Q_day7_tcn_N_aux-T_tc-F_rep1_tech1 n3-d Q day7 tcn N aux-T
n3-d_Q_day7_tcn_N_aux-T_tc-F_rep2_tech1 n3-d Q day7 tcn N aux-T
n3-d_Q_day7_tcn_N_aux-T_tc-F_rep3_tech1 n3-d Q day7 tcn N aux-T
o-d_Q_day7_tcn_N_aux-T_tc-F_rep1_tech1 o-d Q day7 tcn N aux-T
o-d_Q_day7_tcn_N_aux-T_tc-F_rep2_tech1 o-d Q day7 tcn N aux-T
timecourse replicate technical sample_name
<factor> <factor> <factor> <factor>
n3-d_Q_day7_tcn_N_aux-T_tc-F_rep1_tech1 tc-F rep1 tech1 CT8_7716_pIAA_Q_Nascent
n3-d_Q_day7_tcn_N_aux-T_tc-F_rep2_tech1 tc-F rep2 tech1 CT10_7718_pIAA_Q_Nascent
n3-d_Q_day7_tcn_N_aux-T_tc-F_rep3_tech1 tc-F rep3 tech1 CT6_7714_pIAA_Q_Nascent
o-d_Q_day7_tcn_N_aux-T_tc-F_rep1_tech1 tc-F rep1 tech1 CT2_6125_pIAA_Q_Nascent
o-d_Q_day7_tcn_N_aux-T_tc-F_rep2_tech1 tc-F rep2 tech1 CT4_6126_pIAA_Q_Nascent
sizeFactor
<numeric>
n3-d_Q_day7_tcn_N_aux-T_tc-F_rep1_tech1 2.350752
n3-d_Q_day7_tcn_N_aux-T_tc-F_rep2_tech1 1.858782
n3-d_Q_day7_tcn_N_aux-T_tc-F_rep3_tech1 2.518477
o-d_Q_day7_tcn_N_aux-T_tc-F_rep1_tech1 0.333973
o-d_Q_day7_tcn_N_aux-T_tc-F_rep2_tech1 0.289205
# n3-d Q N rep1 1.444025
# n3-d Q N rep2 1.119446
# n3-d Q N rep3 1.547506
# o-d Q N rep1 0.772653
# o-d Q N rep2 0.526895
DESeq2DESeq2 using default parametersdds_SC <- DESeq2::DESeq(dds_SC)
# Check model information
DESeq2::resultsNames(dds_SC)[2]
# [1] "strain_o.d_vs_n3.d" #HERE
#+ Thus, the model varies on strain, OsTIR-depletion is the numerator,
#+ Nab3-depletion is the denominator
#+ - Numerator: "top" in MA plots, "right" in volcano plots
#+ - Denominator: "bottom" in MA plots, "left" in volcano plots
DESeq2::results()Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
Please use `linewidth` instead.
Here, we subset out—i.e., remove—the S. cerevisiae features and are thus only analyzing K. lactis features, using all K. lactis features in the size-factor estimation.
dds_KL <- BiocGenerics::estimateSizeFactors(
dds[dds@rowRanges$genome != "S_cerevisiae", ]
)
dds_KL@colData
# n3-d Q N rep1 0.907373
# n3-d Q N rep2 1.034913
# n3-d Q N rep3 0.944618
# o-d Q N rep1 1.001126
# o-d Q N rep2 1.130992
DESeq2DESeq2 using default parametersdds_KL <- DESeq2::DESeq(dds_KL)
# Check model information
DESeq2::resultsNames(dds_KL)[2]
# [1] "strain_o.d_vs_n3.d" #HERE
#+ Thus, the model varies on strain, OsTIR-depletion is the numerator,
#+ Nab3-depletion is the denominator
#+ - Numerator: "top" in MA plots, "right" in volcano plots
#+ - Denominator: "bottom" in MA plots, "left" in volcano plots
DESeq2::results()
controlGenesTo run analyses using all K.lactis features as
controlGenes, we use a logical vector (a vector composed of
elements with values of either TRUE or FALSE)
obtained from parsing the rowRanges dataframe within the
dds object. In essence, we’re saying, “Return
TRUE if the rowRanges variable
genome has a value of K_lactis; otherwise,
return FALSE.” Then,
BiocGenerics::estimateSizeFactors() is using the values
associated with those TRUEs to isolate the counts for
K. lactis-specific features, in turn using those values to
calculate size factors.
# Recalculate size factors using K. lactis features as `controlGenes`
dds_SC.ctrl_KL <- BiocGenerics::estimateSizeFactors(
dds,
controlGenes = (dds@rowRanges$genome == "K_lactis")
)
dds_SC.ctrl_KL@colData
# Using all of the K. lactis features as `controlGenes`
# n3-d Q N rep1 0.907373
# n3-d Q N rep2 1.034913
# n3-d Q N rep3 0.944618
# o-d Q N rep1 1.001126
# o-d Q N rep2 1.130992
DESeq2DESeq2 using default parametersSince we’ve already calculated the size factors, I think we can exclude K. lactis features from our work from here on out. We have to do some index subsetting to accomplish this (see below).
dds_SC.ctrl_KL <- DESeq2::DESeq(
dds_SC.ctrl_KL[dds_SC.ctrl_KL@rowRanges$genome != "K_lactis", ]
)
# Check model information
DESeq2::resultsNames(dds_SC.ctrl_KL)[2]
# [1] "strain_o.d_vs_n3.d" #HERE
#+ Thus, the model varies on strain, OsTIR-depletion is the numerator,
#+ Nab3-depletion is the denominator
#+ - Numerator: "top" in MA plots, "right" in volcano plots
#+ - Denominator: "bottom" in MA plots, "left" in volcano plots
DESeq2::results()